https://jakubnowosad.com/posts/2024-10-20-spatcomp-bp2/
library(terra)
terra 1.7.23
ndvi2018_tartu = rast("https://github.com/Nowosad/comparing-spatial-patterns-2024/raw/refs/heads/main/data/ndvi2018_tartu.tif")
ndvi2023_tartu = rast("https://github.com/Nowosad/comparing-spatial-patterns-2024/raw/refs/heads/main/data/ndvi2023_tartu.tif")
plot(ndvi2018_tartu, main = "Tartu (2000)")
plot(ndvi2023_tartu, main = "Tartu (2018)")
ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
terra::plot(ndvi_diff)
plot_div = function(r, ...){
r_range = range(values(r), na.rm = TRUE, finite = TRUE)
max_abs = max(abs(r_range))
new_range = c(-max_abs, max_abs)
plot(r, col = hcl.colors(100, palette = "prgn"), range = new_range, ...)
}
plot_div(ndvi_diff)
ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
hist(ndvi_diff)
library(yardstick)
ndvi_rmse = rmse_vec(values(ndvi2023_tartu)[, 1], values(ndvi2018_tartu)[, 1])
ndvi_rmse
[1] 0.2191853
library(diffeR)
Loading required package: ggplot2
ndvi_mad = MAD(ndvi2023_tartu, ndvi2018_tartu)
[1] "The mean of grid1 is less than the mean of grid2"
ndvi_mad$Total
[1] 0.184306
ndvi_cor = focalPairs(c(ndvi2023_tartu, ndvi2018_tartu), w = 5,
fun = "pearson", na.rm = TRUE)
plot_div(ndvi_cor)
library(geodiv)
Attaching package: ‘geodiv’
The following objects are masked from ‘package:terra’:
rotate, sds
window = matrix(1, nrow = 5, ncol = 5)
ndvi2018_tartu_sa_mw = focal_metrics(ndvi2018_tartu, window = window,
metric = "sa", progress = FALSE)
ndvi2023_tartu_sa_mw = focal_metrics(ndvi2023_tartu, window = window,
metric = "sa", progress = FALSE)
ndvi_diff_sa_mw = ndvi2023_tartu_sa_mw$sa - ndvi2018_tartu_sa_mw$sa
plot_div(ndvi_diff_sa_mw)
plot(ndvi2023_tartu)
plot(ndvi2023_tartu_sa_mw$sa)
library(rasterdiv)
ndvi2018_tartu_int = as.int(ndvi2018_tartu * 100)
ndvi2023_tartu_int = as.int(ndvi2023_tartu * 100)
ndvi2018_tartu_rao = paRao(ndvi2018_tartu_int, window = 5, progBar = FALSE)
Warning: Simplify=0. Rounding data to 0 decimal places.
Processing alpha: 1 Moving Window: 5
Processing alpha: 1 Moving Window: 5
Error in utils::combn(p, m = 2, FUN = prod, na.rm = TRUE) : n < m
ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
ndvi_diff_autocor = autocor(ndvi_diff, method = "moran", global = FALSE)
plot_div(ndvi_diff, main = "Difference")
plot_div(ndvi_diff_autocor, main = "Moran's I of the difference")
library(SSIMmap)
ndvi_ssim = ssim_raster(ndvi2018_tartu, ndvi2023_tartu, global = FALSE, w = 3)
plot_div(ndvi_ssim[[1]])
library(waywiser)
cell_sizes = seq(50, 300, by = 50)
ndvi_multi_scale = ww_multi_scale(truth = ndvi2018_tartu, estimate = ndvi2023_tartu,
metrics = list(yardstick::rmse),
cellsize = cell_sizes,
progress = FALSE)
ℹ The package "exactextractr" is required.
✖ Would you like to install it?
1: Yes
2: No
1
Installing package into ‘C:/Users/stahlm/AppData/Local/R/win-library/4.2’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.2/exactextractr_0.10.0.zip'
Content type 'application/zip' length 1448767 bytes (1.4 MB)
downloaded 1.4 MB
package ‘exactextractr’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\stahlm\AppData\Local\Temp\RtmpWCibE9\downloaded_packages
Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.
ndvi_multi_scale
library(diffeR)
ndvi_mad = MAD(ndvi2023_tartu, ndvi2018_tartu, eval = "multiple")
[1] "The mean of grid1 is less than the mean of grid2"
ndvi_mad
library(SSIMmap)
ssim_raster(ndvi2018_tartu, ndvi2023_tartu, global = TRUE)
SSIM: 0.63845 SIM: 0.90432 SIV: 0.90778 SIP: 0.75671